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ABSTRACT 



A method for calculating the flow field in aircraft engine inlets is presented. The 
phenomena of inlet unstart and restart inlet are investigated. Solutions of the 
Reduced Navier Stokes (RNS) equations are obtained with a time consistent direct 
sparse matrix solver that computes the transient flow field both internal and 
external to the inlet. Time varying shocks and time varying recirculation regions 
can be efficiently analyzed. The code is quite general and is suitable for the 
compution of flow for a wide variety of geometries and over a wide range of Mach and 
Reyno Ids numb e r s . 


INTRODUCTION 


At supersonic flight mach numbers the flow within an aircraft engine inlet is 
very complex. Oblique and normal shocks are present and viscous interactions can be 
quite significant. These effects are amplified by the presence of adverse pressure 
gradients, and shock-boundary layer interaction plays an important role in 
determining the flow behaviour. At supersonic flight conditions, the performance of 
the inlet/diffuser depends critically on the nature of the shock interaction and, 'in 
particular, on the location of the terminal shock. At off-design flight conditions, 
for a given diffuser geometry, a shock can move ahead of the cowling so that inlet 
unstart occurs. This causes a sharp reduction in mass flow and pressure recovery and 
an associated large increase in drag. The shock can be swallowed (engine restart) by 
increasing the throat area or by decreasing the back pressure. The performance of 
the inlet then returns to the design value. These complex flow phenomena cannot be 
accurately predicted with analytical techniques. Computational methods can provide a 
reasonable estimate of the critical flow phenomena. Toward this goal, the present 
authors have previously investigated inviscid and low Reynolds numbers flow in two 
dimensional and axisymmetric inlets (1,2,3). 


The present study deals with an investigation of the flow characteristics in 
typical aircraft engine inlets at large Reynolds number. A two dimensional/ 
axisymmetric flow solver for inlet flow field studies has been developed. The 
governing equations are written in general non- orthogonal curvilinear coordinates 
and are discretized using a form of flux vector splitting (4) . This procedure was 
formulated previously by Rubin and Khosla (4,5). The discretization procedure 
remains the same for viscous and inviscid, incompressible, subsonic, transonic and 
supersonic flows, see Rubin and Khosla (6-10), and has been extended for three 
dimensional (11), as well as unsteady flow (1,2,3,12,13) computations. A significant 
feature of this procedure is that it does not require the addition of explicit 
artificial viscosity. Numerical dissipation due to the accuracy of the 
discretization can be minimized on sufficiently fine grids, see Rubin and Himansu 
(14) . 



For unsteady flow computations, in order to capture the transient behaviour 
efficiently, a more robust and time consistent procedure has been developed. A 
direct sparse matrix solver (15,16) has been appropriately modified for coupled 
systems of equations and is applied herein. The choice of the direct solver is 
dictated by considerations of stability, robustness, accuracy and time consistency. 
For steady computations, the solution technique permits large time increments and 
has strong convergence properties; whereas, for transient flows, time consistency is 
maintained in an efficient manner. Implicit, time consistent procedures based on 
approximate factorization, e.g. ADI techniques, typically do not have strong 
convergence properties and may require added transient or steady state artificial 
viscosity (17). The sparse matrix direct solver, considered herein, retains the 
simplicity and robustness of the time marching procedure and, as in the steady-state 
global relaxation formulation, explicitly added artificial viscosity is not 
required. Moreover, for time dependent computations, the time step limitation for 
the direct solver is much less severe than that for other time marching procedures . 


The flow in an inlet at large Reynolds number is inherently turbulent. A simple 
algebraic two layer eddy viscosity (Baldwin- Lomax) model is used herein for 
turbulent flow computations. Section 2 presents the governing equations. The 
boundary conditions and discretization are described in Sections 3 and 4, 
respectively. Section 5 deals with the solution procedure and the results are 
discussed in Section 6. 


GOVERNING EQUATIONS 


The RNS equations are obtained from the full NS equations by neglecting the 
viscous diffusion terms in an appropriate streamwise direction, as well as all 
viscous diffusion terms in the normal momentum and energy balances. In this context, 
the RNS system represents a composite of the Euler plus second-order boundary layer 
equations . The conservation form of the RNS equations are written in general non- 
orthogonal curvilinear co-ordinates (£,r?), so that an arbitrary grid generation 
technique can be applied. 

Continuity Equation: 

(pgr)^ + (pgrU)^ + (pgrV) = 0. (la) 

X Momentum Equation: 

[ Pgr (UX^ + VX )] f + (pgrU 2 X^ + (pgrUVX^ + (pgrUVX^ + (pgrV^)^ + 
r(PY f7 ) e - r(PY^ - (( prx^f U(X^- Y^) + V(XX f Y^)]^ 

-{{ 2/xrY^ [ 2 (UX^+ VX Y^) + (VrY X^ + UrY^) /r] /(3g) y/ (ReJ-0 (lb) 

Y Momentum Equation: 

[pgr(VY f? + Uy] r + (pgrU 2 Y^ + (pgrUVY^ + (pgrUVY^ + (pgrV^)^ 

-r(PX^ + r(PX^ - {{ prY^ U(Y^- X^) + V(Y fj Y ? - X^) ] ^/g) ^ - r{ 2pX^ 

[3(UY^+ VY^X^ - (VrX^ + UrX^yr +(UY^ + VX , Y ^) jj ]/(3 g ) } fj /(Re #o ) - 

(r M Y ? [ U(Y^ - X^) + V(Y^ - /(g) ^ /(Rej=0 (lc) 

Energy Equation: 

(pgrH) r + (pgrHU) ^ + (pgrHV)^ = {r( 7 -l)M 2 / [ 1+ . 5 ( 7 - 1)M 2 ] ) (Pg) ^ - 

{[p rY^HY^/g]^ + [prX^ (HX^) ^/g] ^ )/ (Re^Pr) (Id) 
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Equation of State: 

P + ( 7 -l)pV 2 /(27) - (7-1) [ l/(( 7 -DM 2 )) + .5 ] (pH/ 7 ) (le) 

U and V are the contravariant velocity components in the £ and rj directions, 
respectively; p is the density; P is the pressure; H is the total enthalpy; p is the 
coefficient of viscosity; 7 is the ratio of specific heats; Re^ is the reynolds 

number; Pr is prandtl number and is the free stream mach number. The quantities 


X 


X 


Y* , Y are metrics associated with the coordinate transformation; r is the 

£ n 


V V 

radial location. Two dimensional flow equations are recovered by setting r to one in 
(la-e) ; g is the transformation jacobian and r is time. For viscous flow 
computations , the X and Y momentum equations are appropriately combined to obtain 
the momentum balance in the £ and rj directions. Viscous diffusion terms are retained 
in the £ momentum equation; all viscous dissipation terms in the rj momentum and 
energy equations are neglected. 


In the equations (la-e) , all distances have been normalized with respect to the 
inlet throat radius; the velocities, density, temperature, total enthalpy, viscosity 
are non-dimensionalized with respect to the corresponding freestream values; the 
pressure is non-dimensionalized with respect to twice the free stream dynamic 
pressure . 


INITIAL AND BOUNDARY CONDITIONS 


For viscous flow computations, the inviscid values are prescribed for the 
initial guess. Boundary conditions are such that, at the inflow U, V, p, P and H are 
all prescribed. At the outflow, for boundaries outside of the inlet and far from 
regions of reversed flow, the negative eigenvalue fluxes are neglected. For internal 
flow boundaries, the back pressure is specified. This is consistent with the 
operational or experimental conditions of the inlet. Far from the surface of the 
cowling, uniform flow conditions are imposed. At the surface, for inviscid flow 
calculations, zero normal velocity or injection is specified. For viscous flow 
computations, additional no slip and adiabatic wall temperature conditions have been 
prescribed for most computations. However, cold wall temperature conditions have 
also been considered in some calculations. A wall pressure condition is not 
required. The surface pressure is computed as part of the solution. For external 
outer boundaries, the free stream pressure is specified. More details on boundary 
conditions as applied to internal flows are available in reference (1,2,3,18). 


DISCRETIZATION 


The RNS equations are discretized based on a pressure flux- split technique. The 
differencing has been described in previous references (1-4) and is only briefly 
reviewed here. All convective or £ derivatives are upwind or flux vector 
differenced. The rj derivatives in the continuity and r? momentum equations are mid 
point trapezoidal or two-point (central) differenced, whereas, the rj derivatives in 
the £ momentum and energy equations are three point central differenced. The three 
point t] differencing in the £ momentum and energy equations works quite well for 
normal shocks, but leads to oscillations ahead of strong oblique shocks. These 
oscillations are eliminated with mid-point differencing, similar to that used for 
the continuity and r? momentum equations. For compression regions the £ momentum 
equation, written at an appropriate half point is employed. For expansion regions, 
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the original central differencing scheme is retained. The details of this analysis 
are given in reference (1) . 


The streamwise pressure gradient is flux-split (4). This splitting is consistent 
with the flow physics and does not involve any discontinuous switching across shocks 
or contact discontinuities. The pressure gradient is given by 

P £ = < ' W i-l/2' ) (P i" P i-l' )//A ^i + ^ ' w i+l/2^ P i+l" P i^ A ^i+l 
The parameter w is computed as follows . For unsteady flows , where a differential 
form of the energy equation is employed, see reference (4) , the cartesian form of w 
is 

u> = for < 1 . and w = 1 . for > 1 . 

For curvilinear co-ordinates the eigenvalue analysis indicates that the parameter to 
should be redefined as follows (for details see reference 3) 

w = M? . g^/(Y^ + X^ ) for M. < 1. and u> = 1 for M*. > 1 . 

Where Y , X and g are the metric quantities described previously in section 2. 


The flux form of the streamwise pressure gradient term, with w given at the half 
point is second-order accurate, see reference (4), and captures very sharp normal 
shocks, e.g. over three grid points. It should be noted that the flux splitting is 
employed only in the main flow or £ direction. In the normal, and/or secondary flow 
direction, as described previously, central two or three-point differencing is 
applied. This discretization is capable of capturing very strong normal shocks. 
However, a complete flux-split formulation in both co-ordinate directions. is 
required when considering very strong oblique shock waves, e.g., hypersonic free 
streams. It must be noted that, though the convective streamwise derivatives are 
approximated using first-order differencing the accuracy of the overall scheme is 
somewhere between first and second-order for RNS solutions. This analysis has 
previously been discussed in references (5,8). In reversed flow regions, the 
streamwise convection terms in the £ momentum and energy equation are flux vector or 
upwind differenced and this requires that the parameter w be set to zero, see 
references (2,3,4,11,14). An alternate form of flux- splitting in both £ and r) 
directions is currently being examined. 


In subsonic attached flow regions, upstream influence is transmitted through the 
negative eigenvalue flux or forward differenced part of . At the leading edge, 

upstream influences originate from both the upper and lower surfaces of the cowling. 
This is modelled as an averaged form of the £ momentum equation, written at two half 
points. The details of this procedure are discussed in reference (3). 


Turbulence Model 


For turbulent flow computations the Baldwin-Lomax (B-L) model (19) is used. This is 
an algebraic, two layer eddy viscosity model. The inner model is based on the 
Prandtl-Van Driest formulation and the outer model is based on a modified clauser 
formulation. The distribution of vorticity is used to determine the length scales so 
that the necessity for finding the outer edge of the boundary layer is removed. 
Details of this model are discussed in references (19,20). It is known that for 
highly curved geometries or large regions of recirculation of the B-L model is 
inadequate. For the present study these effects are generally not dominant. 
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SOLUTION PROCEDURE 


The discretized equations are quasilinearized and written in a nine point star 
in delta form. 

E. . 8<j> . ,, . + AM. .86 . , . , + 

ij r i+lj ij r i-lj-l 


A . . 8<j> . . , 
ij iJ-1 


+ B . .6<f>. . 

ij 


+ C . . 8<f> . . 
ij tj+l 


D. . 

ij i-lj 


CM + EM. .6^. , . , +EP..S«L . = G. . 

ij ^l-lj+l ij r i+lj-l ij *i+lj+l ij 

where §</> is the solution vector, the coefficients 

and G. . is a (5x1) vector . 
ij 


( 2 ) 

EM. . are (5x5) matrices 
ij 


For supersonic regions where there is no upstream influence, , EM„ and 

EP. . are zero and the method reduces to a standard initial value problem. This 
ij 

system can then be easily solved using a marching technique such as line relaxation. 
However, for subsonic flow fields and fine meshes the convergence rates of such 
iterative techniques are generally significantly reduced. Moreover, for unsteady 
flows the physical time step limitation can also be severe (13) . Iterative schemes 
for strongly interacting flows are also susceptible to false transients that further 
reduce time step requirements. Most of the difficulties that are associated with 
iterative or approximate factorization techniques can be overcome by the use of a 
direct solver. In the present study the choice of direct solver is dictated by 
considerations of stability, robustness and time consistency. 


The Yale Sparse Matrix Package (YSMP) , developed by Eisenstadt (14) , and 
modified for coupled systems [ 15 ] and f or the boundary conditions detailed 
previously, is applied here. This is an efficient solver as it stores only non-zero 
elements, and reorders the matrix using a minimum degree algorithm to minimize fill- 
in during the LU decomposition. However, for fine meshes and large numbers of mesh 
points, the memory requirement for the direct solver can become significant. 
Although memory associated with present day computers has constantly been 
increasing, access is still limited. To overcome this limitation, a domain 
decomposition strategy is employed. The computational domain is appropriately split 
into subdomains with suitable overlap between adjacent regions. Since we are dealing 
with flow fields that involve moving shocks, the overlap has to be sufficiently 
broad to accomodate the complete shock. Since captured shocks are spread over three 
to four grid points, an overlap of five points is specified. For further details see 
reference (2) . 


RESULTS 


In the first part of the present investigation a two dimensional inlet (flat 
plate diffuser) is considered. The present authors have previously investigated 
inviscid and viscous laminar flow (Re# 12000) for this geometry, see references 
(1,2). Rather interesting flow physics is highlighted even by this simple geometry. 
Flow behaviour within a more complex geometry, axisymmetric inlet with a centerbody 
and cowling was also investigated. The results of this analysis are discussed in 
reference (3). Most results are for = 2.5. Higher Mach numbers were considered 

and these results are presented in reference [2] . 


The unstart and restart of a two dimensional inlet are investigated herein. 
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Supersonic flow between two parallel plates with a sufficiently high back pressure 
is computed. Computations were performed on a 57 x 154 grid (AX= 0.024, AYmin= 
0.00003). In the normal direction, 97 points are placed inside the diffuser and the 
remainder are appropriately distributed outside. The internal flow field (between 
the center-line and cowling) and the external flow field, are computed 
simultaneously. For this grid, the CPU time per iteration on the Cray Y-MP8 is about 
35 seconds per iteration. The memory requirement for the above grid using domain 
decomposition strategy (5 domains each of approximately 15 x 154) is a little over 
three megawords . 


In the present investigation, unstart and restart are investigated by changing 
the back pressure. At unstart a bow shock stands ahead of the inlet (fig.la-b). An 




Fig. lb Unstarted Inlet 


inviscid solution, for an unstarted inlet, as described in reference (1) , was used 
as an initial solution to compute the viscous flow. Steady viscous turbulent flow 
(Re# 1000000) solutions for an unstarted inlet were then obtained. Convergence at 
each time step requires approximately 3 to 5 iterations or 150 CPU seconds. The 
steady state viscous results requires 25 time steps. The unstarted inlet (back 
pressure ratio Pb/Pi=7.2) was restarted by decreasing the back pressure ratio to 6. 
This is well below the critical value. Figures 2a-b depict the time history of the 




Fig. 2a Time History of Unstart 
(Centerline Pressure) 


Fig. 2b Time History of Unstart 
(Cowling Pressure) 
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restart process. Once the shock is swallowed the adverse pressure gradient 
af^jsociated with the shock generates a reversed flow region at the cowling surface. 
S/ince the shock location varies with time, the separation region is also in constant 
flux. A steady state is never achieved. 


tu 

oo 

At a non-dimensionalized time r of 5.0 where r = — — and L is the throat half- 


width, the normal shock is swallowed (inlet restart) and lies very close to the 
inlet tip. At this location the terminal normal shock merges with the inlet leading 
edge oblique shock. The cowling pressure distribution (fig. 2b) indicates a single 
shock close to the inlet lip. However, at time 5.54 the pressure distribution 
(fig. 2b) indicates a double shock. Due to the separation, the boundary layer 
thickens considerably and this generates a strong leading edge oblique shock. At 
time 7.7 the pressure contours (fig. 3a) indicate a very complex shock structure 
close to the surface. An oblique shock originates from the leading edge of the 
cowling lip and this shock intersects the triple point associated with the normal 
shock. The velocity vectors (fig. 3b) depict the recirculation zone associated with 
the triple point shock. 



Fig. 3a Shock Boundary Layer 


Fig. 3b Velocity Vectors 


Interaction (pressure contours) 


Unstart of the inlet was then obtained by gradually increasing the back pressure 
from 6. to 7.2. As a first step, a restarted inlet with the swallowed shock very 
close to the cowling lip was considered. Laminar flow at a Reynolds number of 12000 
was first investigated. Figures 4a-b depict the time history leading to unstart. 




Fig. 4a Time History of Unstart 
(mach# 2.5, Re# 12000) 


Fig. 4b Time History of Unstart 
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Time history is recorded from the time instant (t==0) at which the back pressure is 
held constant at Pb/Pi=7.2. Figures 4c -d indicate that the pressure along the 
cowling at the shock gradually rises and a pressure pulse travels within the 
boundary layer towards the inlet exit. This pressure pulse is associated with a shed 
vortex. The skin friction, figs. 4c-d, and the velocity vector plots (figs. 5a-b) 




Fig. 4c Skin Friction Fig.4d Skin Friction 



Fig. 5a Velocity Vectors 


J.ORO 
1.040 
1.020 
1000 
0.980 
0 900 


Velocity Vectors {unstart lima- 4.0) 
Mach# 2.5. Raff 12000. Pb/PI~7.2 


i "} T r 



Q.3S0 



Fig. 5b Velocity Vectors 


T 


depict the recirculation region associated with the moving vortex. The vortex 
diminishes in size and strength as it moves towards the inlet exit. These results 
indicate that during unstart as the shock moves towards the cowling lip, vorticity 
is shed and then convected downstream with the flow until it is dissipated within 
the boundary- layer . The resulting flow phenomena associated with unstart is 
therefore quite complex as it involves moving shocks and recirculation regions as 
well as vortex shedding and associated viscous -inviscid interactions. Current 
investigations include the effect of surface bleed to reduce the severity of vortex 
interaction. 


The effects of the turbulence model was then considered. Unstart flow conditions 
obtained for a Reynolds number of 1000000 were investigated. The time history 
leading from start to unstart is depicted in Figures 6a-b. A comparison of Figure 2b 
(restart) and Figure 6a (unstart) indicate that the shock speed for the unstart 
condition is much smaller than that for the restart condition. The skin friction 
(fig. 6b) indicates that as the shock moves out of the inlet, the associated 
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Fig. 6a Time History of TJnstart 
(mach# 2.5, Re=1000,000) 


recirculation bubble diminishes in : 
(time=19.4) conditions are achieved, 
vorticity is immediately dissipated and 



Fig. 6b Skin friction 


ize and eventually vanishes when unstart 
For the turbulent flow calculations , the 
shedding is not observed. 


For the flat plate diffuser, which represents a starting point for the 
simulation of transient inlet behaviour, neither computational or experimental 
results are available for comparison. However, there is qualitative agreement with 
flow visualization studies. This has been discussed in reference (2). The only 
detailed and quantitative assessment of accuracy is through grid refinement studies. 
These have previously been carried out and are discussed in reference (2) . The code 
has also been validated for other geometries by comparing with other known solutions 
in reference (1) . 


In the second part of this study, supersonic flow in an axisymmetric inlet with 
a centerbody (fig. 7) is investigated. In this study, the inviscid flow solution has 
been obtained. Computations were performed on an 89 x 115 grid. For supersonic 
inflow, a conical shock is formed at the tip of the centerbody. A second shock 
originates from the tip of the cowling and lies within the inlet (fig. 8) . This 




Fig. 7 Axisymmetric Inlet 


Fig. 8 Pressure Contours 


represents the design case and there is no mass spillage. If the back pressure at 
the inlet exit is raised, a terminal normal shock is formed. Figure 9a shows the 
pressure distribution for a back pressure to free stream pressure ratio (Pb/Pi) of 



8. A terminal shock is formed and lies close to the exit (fig. 9b). If the back 
pressure is sufficiently high, the mass flow behind the normal shock is larger than 
that allowable for the given throat area. As a result, the shock moves toward the 
throat of the inlet and is eventually expelled from the inlet. This results in 
spillage of the excess mass and inlet unstart. This phenomena is further 
investigated herein. 


A sufficiently large back pressure is applied; i.e., (Pb/Pi) = 9. The time 
history leading to unstart is depicted in figures 9a-b. The pressure contours (figs 




Fig. 9a Pressure Distribution 


Fig. 9b Pressure Contours 


10a- d) depict the complex shock pattern that exists at various time intervals during 
unstart. At a non-dimensional time of 208., the terminal normal shock interacts with 
the cowling conical shock (fig. 10c). At a time of 259, the shock is expelled (fig. 
lOd) and the inlet unstarts. The expelled shock interacts with the centerbody 
conical shock and this results in a bow shock ahead of the inlet (fig.lOd). The 
surface pressure distribution at various times is shown in figures lla-b (fig. 10). 
Viscous computations for this geometry are under investigation and these results 
will be reported in a subsequent paper. 




Fig. 10a Pressure Contours Fig. 10b Pressure Contours 
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Fig. 10c Pressure Contours Fig.lOd Pressure Contours 


Centerbody pressure, lnvlscfd, mach#2.S 



Cowling pressure, Invlsdd. Pb/PI-9, Mach#2.5 



Fig. 11a Centerbody Pressure 
Variation 


Fig. lib Cowling Pressure Variation 


SUMMARY 


A RNS flux-split computational procedure has been applied to investigate 
unstart/restart of aircraft engine inlets. A two-dimensional flat plate diffuser and 
an axisymmetric inlet with a centerbody has been considered. A sparse matrix direct 
solver combined with a domain decomposition strategy has been used to efficiently 
compute the complete transient flow behavior for a variety of back pressure ratios. 


Unstart and/or restart of the inlet was initiated by changing the back pressure. 
The associated transients were efficiently captured. The applicability of the RNS 
flux-split procedure for unsteady flows involving moving shocks, time varying 
recirculation regions and shed vorticity has been demonstrated. 
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